MOFA factors

load data, libraries, functions

library(tidyverse)

## Warning: package 'ggplot2' was built under R version 4.3.2

## Warning: package 'tidyr' was built under R version 4.3.2

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(MOFA2)

## Warning: package 'MOFA2' was built under R version 4.3.2

## 
## Attaching package: 'MOFA2'
## 
## The following object is masked from 'package:stats':
## 
##     predict

library(ggplot2)
library(MOFAdata)

## models

load('../data/regression_models/regression_models_IgG_PT_1.RData')
models[[1]]$model_coef

##     Day0.Factor11     Day0.Meta.age Day0.TcellPol.pol       (Intercept) 
##        -0.7229914        -0.2739348        -0.1749157        11.8011841

Day0 T cell pol

Age

Factor 11

MOFAobject = load_model('../data/MOFA_models/MOFA2_noScale_train2020_21.hdf5')
MOFAobject

## Trained MOFA with the following characteristics: 
##  Number of views: 4 
##  Views names: geneExp ab cytokineL cell_freq 
##  Number of features (per view): 6650 31 14 39 
##  Number of groups: 1 
##  Groups names: group1 
##  Number of samples (per group): 508 
##  Number of factors: 15

plot_variance_explained(MOFAobject, x="view", y="factor")

plot_factor(MOFAobject, 
            factor = 11,
            color_by = "IgG_PT"
)

plot_weights(MOFAobject,
             #view = "geneExp",
             factor = 11,
             nfeatures = 20,     # Number of features to highlight
             scale = T,          # Scale weights from -1 to 1
             abs = F             # Take the absolute value?
)

plot_weights(MOFAobject,
             view = "geneExp",
             factor = 11,
             nfeatures = 10,     # Number of features to highlight
             scale = T,          # Scale weights from -1 to 1
             abs = F             # Take the absolute value?
)

plot_weights(MOFAobject,
             view = "cell_freq",
             factor = 11,
             nfeatures = 10,     # Number of features to highlight
             scale = T,          # Scale weights from -1 to 1
             abs = F             # Take the absolute value?
)

plot_weights(MOFAobject,
             view = "ab",
             factor = 11,
             nfeatures = 10,     # Number of features to highlight
             scale = T,          # Scale weights from -1 to 1
             abs = F             # Take the absolute value?
)

plot_weights(MOFAobject,
             view = "cytokineL",
             factor = 11,
             nfeatures = 10,     # Number of features to highlight
             scale = T,          # Scale weights from -1 to 1
             abs = F             # Take the absolute value?
)

plot_top_weights(MOFAobject,
                 #view = "geneExp",
                 factor = 11,
                 nfeatures = 40
)

plot_top_weights(MOFAobject,
                 view = "cell_freq",
                 factor = 11,
                 nfeatures = 40
)

plot_data_heatmap(MOFAobject,
                  view = "geneExp",         # view of interest
                  factor = 11,             # factor of interest
                  features = 60,          # number of features to plot (they are selected by weight)
                  
                  # extra arguments that are passed to the `pheatmap` function
                  cluster_rows = TRUE, cluster_cols = FALSE,
                  show_rownames = TRUE, show_colnames = FALSE
)

plot_data_scatter(MOFAobject,
                  view = "geneExp",         # view of interest
                  factor = 11,             # factor of interest
                  features = 11,           # number of features to plot (they are selected by weight)
                  add_lm = TRUE,         # add linear regression
                  color_by = "Meta.infancy_vac"
)

plot_data_scatter(MOFAobject,
                  view = "cytokineL",         # view of interest
                  factor = 11,             # factor of interest
                  features = 11,           # number of features to plot (they are selected by weight)
                  add_lm = TRUE,         # add linear regression
                  color_by = "Meta.infancy_vac"
)

factors <- get_factors(MOFAobject, factors = "all")
lapply(factors,dim)

## $group1
## [1] 508  15

weights <- get_weights(MOFAobject, views = "all", factors = "Factor11")
a= weights$geneExp
ggplot(a, aes(x=Factor11)) +
    geom_density() 

a['Expr.ENSG00000277632.1', ]

## [1] 0.1655906

a = a[order(a[,1]), ]
head(a)

## Expr.ENSG00000156642.16 Expr.ENSG00000134970.13  Expr.ENSG00000113070.7 
##              -0.3186934              -0.3022599              -0.2715794 
##  Expr.ENSG00000146281.5  Expr.ENSG00000178726.6 Expr.ENSG00000021300.13 
##              -0.2678805              -0.2573302              -0.2512104

tail(a)

##  Expr.ENSG00000213190.3 Expr.ENSG00000092148.12 Expr.ENSG00000006459.10 
##               0.2711598               0.2753295               0.2829757 
## Expr.ENSG00000167548.14 Expr.ENSG00000026652.13 Expr.ENSG00000100403.11 
##               0.2889911               0.3086037               0.3112866

tail(names(a), n=200)

##   [1] "Expr.ENSG00000185000.11" "Expr.ENSG00000197119.12"
##   [3] "Expr.ENSG00000152818.18" "Expr.ENSG00000173889.15"
##   [5] "Expr.ENSG00000120868.13" "Expr.ENSG00000066777.8" 
##   [7] "Expr.ENSG00000112182.14" "Expr.ENSG00000005249.12"
##   [9] "Expr.ENSG00000185920.15" "Expr.ENSG00000065883.14"
##  [11] "Expr.ENSG00000165959.11" "Expr.ENSG00000109906.13"
##  [13] "Expr.ENSG00000185129.5"  "Expr.ENSG00000196923.13"
##  [15] "Expr.ENSG00000160326.13" "Expr.ENSG00000078177.13"
##  [17] "Expr.ENSG00000196313.11" "Expr.ENSG00000145050.15"
##  [19] "Expr.ENSG00000135362.13" "Expr.ENSG00000214021.15"
##  [21] "Expr.ENSG00000007202.14" "Expr.ENSG00000166164.15"
##  [23] "Expr.ENSG00000171467.15" "Expr.ENSG00000197102.10"
##  [25] "Expr.ENSG00000173575.20" "Expr.ENSG00000248333.8" 
##  [27] "Expr.ENSG00000176953.12" "Expr.ENSG00000198399.14"
##  [29] "Expr.ENSG00000221968.8"  "Expr.ENSG00000197976.11"
##  [31] "Expr.ENSG00000165801.9"  "Expr.ENSG00000124782.19"
##  [33] "Expr.ENSG00000277632.1"  "Expr.ENSG00000158545.15"
##  [35] "Expr.ENSG00000073350.13" "Expr.ENSG00000143013.12"
##  [37] "Expr.ENSG00000197535.14" "Expr.ENSG00000100097.11"
##  [39] "Expr.ENSG00000198223.16" "Expr.ENSG00000005339.14"
##  [41] "Expr.ENSG00000185215.8"  "Expr.ENSG00000137073.21"
##  [43] "Expr.ENSG00000076770.14" "Expr.ENSG00000007237.18"
##  [45] "Expr.ENSG00000108175.16" "Expr.ENSG00000154359.12"
##  [47] "Expr.ENSG00000114933.15" "Expr.ENSG00000171604.11"
##  [49] "Expr.ENSG00000141564.13" "Expr.ENSG00000153317.14"
##  [51] "Expr.ENSG00000167258.13" "Expr.ENSG00000145246.13"
##  [53] "Expr.ENSG00000185340.15" "Expr.ENSG00000266173.6" 
##  [55] "Expr.ENSG00000272391.5"  "Expr.ENSG00000134982.16"
##  [57] "Expr.ENSG00000113719.15" "Expr.ENSG00000175787.16"
##  [59] "Expr.ENSG00000109381.19" "Expr.ENSG00000180530.10"
##  [61] "Expr.ENSG00000047365.11" "Expr.ENSG00000137331.11"
##  [63] "Expr.ENSG00000181788.3"  "Expr.ENSG00000120071.13"
##  [65] "Expr.ENSG00000132465.10" "Expr.ENSG00000157106.16"
##  [67] "Expr.ENSG00000184863.10" "Expr.ENSG00000179348.11"
##  [69] "Expr.ENSG00000261609.5"  "Expr.ENSG00000005206.16"
##  [71] "Expr.ENSG00000092439.13" "Expr.ENSG00000067365.14"
##  [73] "Expr.ENSG00000196912.12" "Expr.ENSG00000166860.2" 
##  [75] "Expr.ENSG00000105738.10" "Expr.ENSG00000100307.12"
##  [77] "Expr.ENSG00000146232.15" "Expr.ENSG00000162298.18"
##  [79] "Expr.ENSG00000186814.13" "Expr.ENSG00000055609.17"
##  [81] "Expr.ENSG00000099991.17" "Expr.ENSG00000042980.12"
##  [83] "Expr.ENSG00000140287.10" "Expr.ENSG00000145990.10"
##  [85] "Expr.ENSG00000169020.9"  "Expr.ENSG00000083223.17"
##  [87] "Expr.ENSG00000100888.12" "Expr.ENSG00000232810.3" 
##  [89] "Expr.ENSG00000131127.13" "Expr.ENSG00000161835.10"
##  [91] "Expr.ENSG00000263264.1"  "Expr.ENSG00000116852.14"
##  [93] "Expr.ENSG00000100321.14" "Expr.ENSG00000127603.24"
##  [95] "Expr.ENSG00000182362.13" "Expr.ENSG00000152127.8" 
##  [97] "Expr.ENSG00000119242.8"  "Expr.ENSG00000165671.19"
##  [99] "Expr.ENSG00000116649.9"  "Expr.ENSG00000168159.11"
## [101] "Expr.ENSG00000132680.10" "Expr.ENSG00000167978.16"
## [103] "Expr.ENSG00000173200.12" "Expr.ENSG00000060237.16"
## [105] "Expr.ENSG00000174501.14" "Expr.ENSG00000155657.26"
## [107] "Expr.ENSG00000133812.15" "Expr.ENSG00000144840.8" 
## [109] "Expr.ENSG00000116539.12" "Expr.ENSG00000100485.11"
## [111] "Expr.ENSG00000119778.14" "Expr.ENSG00000090924.14"
## [113] "Expr.ENSG00000174373.16" "Expr.ENSG00000185619.18"
## [115] "Expr.ENSG00000065613.13" "Expr.ENSG00000148187.17"
## [117] "Expr.ENSG00000135913.10" "Expr.ENSG00000160179.18"
## [119] "Expr.ENSG00000170476.15" "Expr.ENSG00000145088.8" 
## [121] "Expr.ENSG00000104517.12" "Expr.ENSG00000162804.13"
## [123] "Expr.ENSG00000106948.16" "Expr.ENSG00000162402.13"
## [125] "Expr.ENSG00000243335.9"  "Expr.ENSG00000169718.17"
## [127] "Expr.ENSG00000090061.17" "Expr.ENSG00000275302.1" 
## [129] "Expr.ENSG00000140090.17" "Expr.ENSG00000004700.15"
## [131] "Expr.ENSG00000039523.19" "Expr.ENSG00000137040.9" 
## [133] "Expr.ENSG00000124203.6"  "Expr.ENSG00000204397.7" 
## [135] "Expr.ENSG00000160305.17" "Expr.ENSG00000137992.14"
## [137] "Expr.ENSG00000011454.16" "Expr.ENSG00000116191.17"
## [139] "Expr.ENSG00000221869.4"  "Expr.ENSG00000135976.17"
## [141] "Expr.ENSG00000100201.20" "Expr.ENSG00000084112.14"
## [143] "Expr.ENSG00000011258.15" "Expr.ENSG00000175727.13"
## [145] "Expr.ENSG00000134215.15" "Expr.ENSG00000128815.19"
## [147] "Expr.ENSG00000184232.8"  "Expr.ENSG00000183530.13"
## [149] "Expr.ENSG00000169385.2"  "Expr.ENSG00000185278.15"
## [151] "Expr.ENSG00000175857.8"  "Expr.ENSG00000160584.15"
## [153] "Expr.ENSG00000104093.13" "Expr.ENSG00000214029.4" 
## [155] "Expr.ENSG00000182325.10" "Expr.ENSG00000160218.12"
## [157] "Expr.ENSG00000184897.5"  "Expr.ENSG00000168769.13"
## [159] "Expr.ENSG00000135837.15" "Expr.ENSG00000149485.18"
## [161] "Expr.ENSG00000182621.17" "Expr.ENSG00000182095.14"
## [163] "Expr.ENSG00000123933.16" "Expr.ENSG00000169252.5" 
## [165] "Expr.ENSG00000173064.12" "Expr.ENSG00000070476.14"
## [167] "Expr.ENSG00000145819.15" "Expr.ENSG00000204220.11"
## [169] "Expr.ENSG00000095564.13" "Expr.ENSG00000134545.13"
## [171] "Expr.ENSG00000130856.15" "Expr.ENSG00000196449.3" 
## [173] "Expr.ENSG00000144711.13" "Expr.ENSG00000198563.13"
## [175] "Expr.ENSG00000099331.13" "Expr.ENSG00000132017.10"
## [177] "Expr.ENSG00000162783.10" "Expr.ENSG00000154237.12"
## [179] "Expr.ENSG00000162775.14" "Expr.ENSG00000177463.15"
## [181] "Expr.ENSG00000153898.12" "Expr.ENSG00000131149.18"
## [183] "Expr.ENSG00000184207.8"  "Expr.ENSG00000135723.13"
## [185] "Expr.ENSG00000188822.7"  "Expr.ENSG00000178951.8" 
## [187] "Expr.ENSG00000198198.15" "Expr.ENSG00000160299.16"
## [189] "Expr.ENSG00000170871.11" "Expr.ENSG00000187837.3" 
## [191] "Expr.ENSG00000082805.19" "Expr.ENSG00000135363.11"
## [193] "Expr.ENSG00000158457.5"  "Expr.ENSG00000198728.10"
## [195] "Expr.ENSG00000213190.3"  "Expr.ENSG00000092148.12"
## [197] "Expr.ENSG00000006459.10" "Expr.ENSG00000167548.14"
## [199] "Expr.ENSG00000026652.13" "Expr.ENSG00000100403.11"

b= weights$ab
ggplot(b, aes(x=Factor11)) +
    geom_density() 

c= weights$cytokineL
ggplot(c, aes(x=Factor11)) +
    geom_density() 

d= weights$cell_freq
ggplot(d, aes(x=Factor11)) +
    geom_density() 

enrichment

against Hallmark

MOFAobject3 = MOFAobject
ens_map = read_tsv('~/Documents/Projects/CMI-PB2/CMI-PB-Oct2023-FINAL/data/raw_prediction_dataset/gene_90_38_export.tsv')

## Rows: 63971 Columns: 12
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (6): biotype, versioned_ensembl_gene_id, description, display_label, syn...
## dbl (6): seq_region_strand, seq_region_start, seq_region_end, canonical_tran...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(ens_map)

## # A tibble: 6 × 12
##   biotype        seq_region_strand seq_region_start seq_region_end
##   <chr>                      <dbl>            <dbl>          <dbl>
## 1 miRNA                         -1         55372940       55373034
## 2 snRNA                         -1         59450673       59450772
## 3 miRNA                          1         41691585       41691678
## 4 protein_coding                 1         54201473       54208260
## 5 misc_RNA                       1         10287413       10287482
## 6 snRNA                          1        109992325      109992431
## # ℹ 8 more variables: versioned_ensembl_gene_id <chr>, description <chr>,
## #   canonical_transcript_id <dbl>, display_label <chr>, synonym <chr>,
## #   dbprimary_acc <dbl>, gene_id <dbl>, gene_summary <chr>

genes = as.character(rownames(MOFAobject3@data$geneExp$group1))
genes = gsub('Expr.', '', genes)
table(genes %in% ens_map$versioned_ensembl_gene_id)

## 
## TRUE 
## 6650

ens_map = ens_map[!duplicated(ens_map$versioned_ensembl_gene_id), ]
rownames(ens_map) = ens_map$versioned_ensembl_gene_id

## Warning: Setting row names on a tibble is deprecated.

genes_sym = ens_map[genes, 'display_label']$display_label

rownames(MOFAobject3@data$geneExp$group1) = genes_sym
names(MOFAobject3@intercepts$geneExp$group1) = genes_sym

MOFAobject3@features_metadata$feature[1:6650] = genes_sym

rownames(MOFAobject3@expectations$W$geneExp) = genes_sym

hallmark = fgsea::gmtPathways('~/Documents/References/Genesets/MSigDB/h.all.v2023.2.Hs.symbols.gmt.txt')


# Get all unique genes
all_genes <- sort(unique(unlist(hallmark)))

# Create the binary matrix
binary_matrix <- sapply(all_genes, function(g) {
    as.numeric(sapply(hallmark, function(x) g %in% x))
})

# Transpose and convert to a data frame for better visualization
#binary_matrix <- as.data.frame(t(binary_matrix))
binary_matrix[1:5, 1:5]

##      A2M AAAS AADAT AARS1 ABAT
## [1,]   0    0     0     0    0
## [2,]   0    0     0     1    0
## [3,]   0    0     0     0    0
## [4,]   0    0     0     0    0
## [5,]   0    0     0     0    0

# Add row names
rownames(binary_matrix) <- names(hallmark)

# Display the binary matrix
binary_matrix[1:5, 1:10]

##                              A2M AAAS AADAT AARS1 ABAT ABCA1 ABCA2 ABCA3 ABCA4
## HALLMARK_ADIPOGENESIS          0    0     0     0    0     1     0     0     0
## HALLMARK_ALLOGRAFT_REJECTION   0    0     0     1    0     0     0     0     0
## HALLMARK_ANDROGEN_RESPONSE     0    0     0     0    0     0     0     0     0
## HALLMARK_ANGIOGENESIS          0    0     0     0    0     0     0     0     0
## HALLMARK_APICAL_JUNCTION       0    0     0     0    0     0     0     0     0
##                              ABCA5
## HALLMARK_ADIPOGENESIS            0
## HALLMARK_ALLOGRAFT_REJECTION     0
## HALLMARK_ANDROGEN_RESPONSE       0
## HALLMARK_ANGIOGENESIS            0
## HALLMARK_APICAL_JUNCTION         0

dim(binary_matrix)

## [1]   50 4384

enrichment3.parametric.pos <- run_enrichment(MOFAobject3,
                                            view = "geneExp", factors = 11,
                                            feature.sets = as.matrix(binary_matrix),
                                            sign = "positive",
                                            statistical.test = "parametric"
)

## Intersecting features names in the model and the gene set annotation results in a total of 2153 features.

## 
## Running feature set Enrichment Analysis with the following options...
## View: geneExp 
## Number of feature sets: 48 
## Set statistic: mean.diff 
## Statistical test: parametric

## Subsetting weights with positive sign

## 

plot_enrichment(enrichment3.parametric.pos, 
                factor = 'Factor11', 
                max.pathways = 15
)

plot_enrichment_detailed(enrichment3.parametric.pos, 
                         factor = 'Factor11', 
                         max.genes = 8, 
                         max.pathways = 5
)

## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

enrichment3.parametric.neg <- run_enrichment(MOFAobject3,
                                             view = "geneExp", factors = 11,
                                             feature.sets = as.matrix(binary_matrix),
                                             sign = "negative",
                                             statistical.test = "parametric"
)

## Intersecting features names in the model and the gene set annotation results in a total of 2153 features.

## 
## Running feature set Enrichment Analysis with the following options...
## View: geneExp 
## Number of feature sets: 48 
## Set statistic: mean.diff 
## Statistical test: parametric

## Subsetting weights with negative sign

## 

against BTM

btm = fgsea::gmtPathways('~/Documents/References/Genesets/BTM/BTM_for_GSEA_20131008.gmt')
lapply(btm, length)

## $`targets of FOSL1/2 (M0)`
## [1] 12
## 
## $`integrin cell surface interactions (I) (M1.0)`
## [1] 29
## 
## $`integrin cell surface interactions (II) (M1.1)`
## [1] 12
## 
## $`extracellular matrix (I) (M2.0)`
## [1] 30
## 
## $`extracellular matrix (II) (M2.1)`
## [1] 45
## 
## $`extracellular matrix (III) (M2.2)`
## [1] 14
## 
## $`regulation of signal transduction (M3)`
## [1] 47
## 
## $`cell cycle and transcription (M4.0)`
## [1] 335
## 
## $`cell cycle (I) (M4.1)`
## [1] 145
## 
## $`PLK1 signaling events (M4.2)`
## [1] 34
## 
## $`myeloid cell enriched receptors and transporters (M4.3)`
## [1] 31
## 
## $`mitotic cell cycle - DNA replication (M4.4)`
## [1] 30
## 
## $`mitotic cell cycle in stimulated CD4 T cells (M4.5)`
## [1] 35
## 
## $`cell division in stimulated CD4 T cells (M4.6)`
## [1] 20
## 
## $`mitotic cell cycle (M4.7)`
## [1] 21
## 
## $`cell division - E2F transcription network (M4.8)`
## [1] 19
## 
## $`mitotic cell cycle in stimulated CD4 T cells (M4.9)`
## [1] 16
## 
## $`cell cycle (II) (M4.10)`
## [1] 14
## 
## $`mitotic cell cycle in stimulated CD4 T cells (M4.11)`
## [1] 12
## 
## $`C-MYC transcriptional network (M4.12)`
## [1] 12
## 
## $`cell junction (GO) (M4.13)`
## [1] 11
## 
## $`Rho GTPase cycle (M4.14)`
## [1] 10
## 
## $`enriched in monocytes (I) (M4.15)`
## [1] 11
## 
## $`regulation of antigen presentation and immune response (M5.0)`
## [1] 81
## 
## $`T cell activation and signaling (M5.1)`
## [1] 25
## 
## $`mitotic cell division (M6)`
## [1] 32
## 
## $`enriched in T cells (I) (M7.0)`
## [1] 62
## 
## $`T cell activation (I) (M7.1)`
## [1] 50
## 
## $`enriched in NK cells (I) (M7.2)`
## [1] 49
## 
## $`T cell activation (II) (M7.3)`
## [1] 31
## 
## $`T cell activation (III) (M7.4)`
## [1] 14
## 
## $`E2F transcription factor network (M8)`
## [1] 18
## 
## $`B cell development (M9)`
## [1] 11
## 
## $`E2F1 targets (Q3) (M10.0)`
## [1] 33
## 
## $`E2F1 targets (Q4) (M10.1)`
## [1] 21
## 
## $`enriched in monocytes (II) (M11.0)`
## [1] 189
## 
## $`blood coagulation (M11.1)`
## [1] 22
## 
## $`formyl peptide receptor mediated neutrophil response (M11.2)`
## [1] 10
## 
## $`CD28 costimulation (M12)`
## [1] 10
## 
## $`innate activation by cytosolic DNA sensing (M13)`
## [1] 11
## 
## $`T cell differentiation (M14)`
## [1] 12
## 
## $`Ran mediated mitosis (M15)`
## [1] 13
## 
## $`TLR and inflammatory signaling (M16)`
## [1] 45
## 
## $`Hox cluster I (M17.0)`
## [1] 16
## 
## $`Hox cluster II (M17.1)`
## [1] 12
## 
## $`Hox cluster III (M17.2)`
## [1] 12
## 
## $`Hox cluster IV (M17.3)`
## [1] 10
## 
## $`T cell differentiation via ITK and PKC (M18)`
## [1] 11
## 
## $`T cell differentiation (Th2) (M19)`
## [1] 17
## 
## $`AP-1 transcription factor network (M20)`
## [1] 15
## 
## $`cell adhesion (lymphocyte homing) (M21)`
## [1] 10
## 
## $`mismatch repair (I) (M22.0)`
## [1] 27
## 
## $`mismatch repair (II) (M22.1)`
## [1] 13
## 
## $`RA, WNT, CSF receptors network (monocyte) (M23)`
## [1] 12
## 
## $`cell activation (IL15, IL23, TNF) (M24)`
## [1] 15
## 
## $`TLR8-BAFF network (M25)`
## [1] 10
## 
## $`TBA (M26.0)`
## [1] 30
## 
## $`TBA (M26.1)`
## [1] 22
## 
## $`TBA (M26.2)`
## [1] 14
## 
## $`chemokine cluster (I) (M27.0)`
## [1] 26
## 
## $`chemokine cluster (II) (M27.1)`
## [1] 15
## 
## $`antigen presentation (lipids and proteins) (M28)`
## [1] 11
## 
## $`proinflammatory cytokines and chemokines (M29)`
## [1] 10
## 
## $`cell movement, Adhesion & Platelet activation (M30)`
## [1] 20
## 
## $`cell cycle and growth arrest (M31)`
## [1] 12
## 
## $`platelet activation (I) (M32.0)`
## [1] 23
## 
## $`platelet activation (II) (M32.1)`
## [1] 21
## 
## $`CORO1A-DEF6 network (I) (M32.2)`
## [1] 20
## 
## $`KLF12 targets network (M32.3)`
## [1] 17
## 
## $`CORO1A-DEF6 network (II) (M32.4)`
## [1] 15
## 
## $`TBA (M32.5)`
## [1] 15
## 
## $`TBA (M32.6)`
## [1] 13
## 
## $`TBA (M32.7)`
## [1] 11
## 
## $`cytoskeletal remodeling (M32.8)`
## [1] 10
## 
## $`inflammatory response (M33)`
## [1] 11
## 
## $`cytoskeletal remodeling (enriched for SRF targets) (M34)`
## [1] 10
## 
## $`signaling in T cells (I) (M35.0)`
## [1] 14
## 
## $`signaling in T cells (II) (M35.1)`
## [1] 11
## 
## $`T cell surface, activation (M36)`
## [1] 12
## 
## $`immune activation - generic cluster (M37.0)`
## [1] 347
## 
## $`enriched in neutrophils (I) (M37.1)`
## [1] 52
## 
## $`endoplasmic reticulum (M37.2)`
## [1] 19
## 
## $`cell division (M37.3)`
## [1] 10
## 
## $`chemokines and receptors (M38)`
## [1] 11
## 
## $`integrin mediated leukocyte migration (M39)`
## [1] 10
## 
## $`complement and other receptors in DCs (M40)`
## [1] 13
## 
## $`TBA (M41.0)`
## [1] 17
## 
## $`TBA (M41.1)`
## [1] 15
## 
## $`TBA (M41.2)`
## [1] 14
## 
## $`TBA (M41.3)`
## [1] 10
## 
## $`ATF targets network (M41.4)`
## [1] 10
## 
## $`platelet activation (III) (M42)`
## [1] 10
## 
## $`myeloid, dendritic cell activation via NFkB (I) (M43.0)`
## [1] 15
## 
## $`myeloid, dendritic cell activation via NFkB (II) (M43.1)`
## [1] 14
## 
## $`T cell signaling and costimulation (M44)`
## [1] 10
## 
## $`leukocyte activation and migration (M45)`
## [1] 11
## 
## $`cell division (stimulated CD4+ T cells) (M46)`
## [1] 28
## 
## $`enriched in B cells (I) (M47.0)`
## [1] 53
## 
## $`enriched in B cells (II) (M47.1)`
## [1] 43
## 
## $`enriched in B cells (III) (M47.2)`
## [1] 37
## 
## $`enriched in B cells (IV) (M47.3)`
## [1] 16
## 
## $`enriched in B cells (V) (M47.4)`
## [1] 11
## 
## $`TBA (M48)`
## [1] 13
## 
## $`transcription regulation in cell development (M49)`
## [1] 47
## 
## $`CD1 and other DC receptors (M50)`
## [1] 10
## 
## $`cell adhesion (M51)`
## [1] 38
## 
## $`T cell activation (IV) (M52)`
## [1] 13
## 
## $`inflammasome receptors and signaling (M53)`
## [1] 12
## 
## $`BCR signaling (M54)`
## [1] 12
## 
## $`TBA (M55)`
## [1] 12
## 
## $`suppression of MAPK signaling (M56)`
## [1] 12
## 
## $`immuregulation - monocytes, T and B cells (M57)`
## [1] 13
## 
## $`B cell development/activation (M58)`
## [1] 11
## 
## $`CCR1, 7 and cell signaling (M59)`
## [1] 10
## 
## $`lymphocyte generic cluster (M60)`
## [1] 17
## 
## $`enriched in NK cells (II) (M61.0)`
## [1] 24
## 
## $`enriched in NK cells (KIR cluster) (M61.1)`
## [1] 13
## 
## $`enriched in NK cells (receptor activation) (M61.2)`
## [1] 16
## 
## $`T & B cell development, activation (M62.0)`
## [1] 44
## 
## $`enriched for unknown TF motif CTCNANGTGNY (M62.1)`
## [1] 12
## 
## $`regulation of localization (GO) (M63)`
## [1] 12
## 
## $`enriched in activated dendritic cells/monocytes (M64)`
## [1] 17
## 
## $`IL2, IL7, TCR network (M65)`
## [1] 12
## 
## $`TBA (M66)`
## [1] 17
## 
## $`activated dendritic cells (M67)`
## [1] 12
## 
## $`RIG-1 like receptor signaling (M68)`
## [1] 10
## 
## $`enriched in B cells (VI) (M69)`
## [1] 20
## 
## $`TBA (M70.0)`
## [1] 20
## 
## $`TBA (M70.1)`
## [1] 10
## 
## $`enriched in antigen presentation (I) (M71)`
## [1] 18
## 
## $`TBA (M72.0)`
## [1] 24
## 
## $`TBA (M72.1)`
## [1] 17
## 
## $`TBA (M72.2)`
## [1] 12
## 
## $`enriched in monocytes (III) (M73)`
## [1] 12
## 
## $`transcriptional targets of glucocorticoid receptor (M74)`
## [1] 12
## 
## $`antiviral IFN signature (M75)`
## [1] 22
## 
## $`DNA repair (M76)`
## [1] 22
## 
## $`collagen, TGFB family et al (M77)`
## [1] 35
## 
## $`myeloid cell cytokines, metallopeptidases and laminins (M78)`
## [1] 11
## 
## $`TBA (M79)`
## [1] 10
## 
## $`TBA (M80)`
## [1] 20
## 
## $`enriched in myeloid cells and monocytes (M81)`
## [1] 35
## 
## $`signal transduction, plasma membrane (M82)`
## [1] 13
## 
## $`enriched in naive and memory B cells (M83)`
## [1] 10
## 
## $`integrins and cell adhesion (M84)`
## [1] 10
## 
## $`platelet activation and degranulation (M85)`
## [1] 12
## 
## $`chemokines and inflammatory molecules in myeloid cells (M86.0)`
## [1] 19
## 
## $`proinflammatory dendritic cell, myeloid cell response (M86.1)`
## [1] 13
## 
## $`transmembrane transport (I) (M87)`
## [1] 24
## 
## $`leukocyte migration (M88.0)`
## [1] 51
## 
## $`enriched in hepatocyte nuclear factors (I) (M88.1)`
## [1] 13
## 
## $`enriched in hepatocyte nuclear factors (II) (M88.2)`
## [1] 12
## 
## $`putative targets of PAX3 (M89.0)`
## [1] 16
## 
## $`putative targets of PAX3 (M89.1)`
## [1] 10
## 
## $`TBA (M90)`
## [1] 12
## 
## $`adhesion and migration, chemotaxis (M91)`
## [1] 13
## 
## $`lipid metabolism, endoplasmic reticulum (M92)`
## [1] 10
## 
## $`TBA (M93)`
## [1] 10
## 
## $`growth factor induced, enriched in nuclear receptor subfamily 4 (M94)`
## [1] 12
## 
## $`enriched in antigen presentation (II) (M95.0)`
## [1] 24
## 
## $`enriched in antigen presentation (III) (M95.1)`
## [1] 15
## 
## $`Hox cluster V (M96)`
## [1] 10
## 
## $`enriched for SMAD2/3 signaling (M97)`
## [1] 10
## 
## $`TBA (M98.0)`
## [1] 18
## 
## $`TBA (M98.1)`
## [1] 10
## 
## $`TBA (M99)`
## [1] 10
## 
## $`MAPK, RAS signaling (M100)`
## [1] 10
## 
## $`phosphatidylinositol signaling system (M101)`
## [1] 14
## 
## $`TBA (M102)`
## [1] 12
## 
## $`cell cycle (III) (M103)`
## [1] 51
## 
## $`TBA (M104)`
## [1] 10
## 
## $`TBA (M105)`
## [1] 18
## 
## $`nuclear pore complex (M106.0)`
## [1] 17
## 
## $`nuclear pore complex (mitosis) (M106.1)`
## [1] 11
## 
## $`Hox cluster VI (M107)`
## [1] 11
## 
## $`TBA (M108)`
## [1] 11
## 
## $`receptors, cell migration (M109)`
## [1] 15
## 
## $`axon guidance (M110)`
## [1] 10
## 
## $`viral sensing & immunity; IRF2 targets network (I) (M111.0)`
## [1] 17
## 
## $`viral sensing & immunity; IRF2 targets network (II) (M111.1)`
## [1] 11
## 
## $`complement activation (I) (M112.0)`
## [1] 17
## 
## $`complement activation (II) (M112.1)`
## [1] 10
## 
## $`golgi membrane (I) (M113)`
## [1] 10
## 
## $`TBA (M114.0)`
## [1] 39
## 
## $`glycerophospholipid metabolism (M114.1)`
## [1] 11
## 
## $`cytokines - recepters cluster (M115)`
## [1] 11
## 
## $`TBA (M116)`
## [1] 13
## 
## $`cell adhesion (GO) (M117)`
## [1] 21
## 
## $`enriched in monocytes (IV) (M118.0)`
## [1] 57
## 
## $`enriched in monocytes (surface) (M118.1)`
## [1] 17
## 
## $`enriched in activated dendritic cells (I) (M119)`
## [1] 11
## 
## $`TBA (M120)`
## [1] 11
## 
## $`TBA (M121)`
## [1] 12
## 
## $`enriched for cell migration (M122)`
## [1] 11
## 
## $`enriched in B cell differentiation (M123)`
## [1] 13
## 
## $`enriched in membrane proteins (M124)`
## [1] 19
## 
## $`TBA (M125)`
## [1] 11
## 
## $`double positive thymocytes (M126)`
## [1] 13
## 
## $`type I interferon response (M127)`
## [1] 12
## 
## $`TBA (M128)`
## [1] 12
## 
## $`inositol phosphate metabolism (M129)`
## [1] 13
## 
## $`enriched in G-protein coupled receptors (M130)`
## [1] 10
## 
## $`TBA (M131)`
## [1] 14
## 
## $`recruitment of neutrophils (M132)`
## [1] 11
## 
## $`cell adhesion, membrane (M133.0)`
## [1] 14
## 
## $`cell cell adhesion (M133.1)`
## [1] 11
## 
## $`Membrane, ER proteins (M134)`
## [1] 17
## 
## $`enriched in plasma membrane proteins (I) (M135.0)`
## [1] 16
## 
## $`enriched in plasma membrane proteins (II) (M135.1)`
## [1] 15
## 
## $`TBA (M136)`
## [1] 17
## 
## $`TBA (M137)`
## [1] 16
## 
## $`enriched for ubiquitination (M138)`
## [1] 11
## 
## $`lysosomal/endosomal proteins (M139)`
## [1] 11
## 
## $`extracellular matrix, complement (M140)`
## [1] 14
## 
## $`TBA (M141)`
## [1] 27
## 
## $`transmembrane and ion transporters (I) (M142)`
## [1] 13
## 
## $`nuclear pore, transport; mRNA splicing, processing (M143)`
## [1] 11
## 
## $`cell cycle, ATP binding (M144)`
## [1] 16
## 
## $`cytoskeleton/actin (SRF transcription targets) (M145.0)`
## [1] 16
## 
## $`cytoskeleton/actin (SRF transcription targets) (M145.1)`
## [1] 14
## 
## $`MHC-TLR7-TLR8 cluster (M146)`
## [1] 17
## 
## $`intracellular transport (M147)`
## [1] 17
## 
## $`TBA (M148)`
## [1] 11
## 
## $`TBA (M149)`
## [1] 10
## 
## $`innate antiviral response (M150)`
## [1] 12
## 
## $`TBA (M151)`
## [1] 10
## 
## $`TBA (source: B cells) (M152.0)`
## [1] 25
## 
## $`TBA (source: naive B cells) (M152.1)`
## [1] 21
## 
## $`TBA (source: memory B cells) (M152.2)`
## [1] 18
## 
## $`TBA (M153)`
## [1] 16
## 
## $`amino acid metabolishm and transport (M154.0)`
## [1] 33
## 
## $`transmembrane transport (SLC cluster) (M154.1)`
## [1] 17
## 
## $`G protein coupled receptors cluster (M155)`
## [1] 10
## 
## $`plasma cells & B cells, immunoglobulins (M156.0)`
## [1] 43
## 
## $`plasma cells, immunoglobulins (M156.1)`
## [1] 36
## 
## $`enriched in NK cells (III) (M157)`
## [1] 14
## 
## $`interferon alpha response (I) (M158.0)`
## [1] 16
## 
## $`interferon alpha response (II) (M158.1)`
## [1] 13
## 
## $`G protein mediated calcium signaling (M159)`
## [1] 10
## 
## $`leukocyte differentiation (M160)`
## [1] 16
## 
## $`TBA (M161)`
## [1] 14
## 
## $`plasma membrane, cell junction (M162.0)`
## [1] 18
## 
## $`cell junction (M162.1)`
## [1] 11
## 
## $`enriched in neutrophils (II) (M163)`
## [1] 14
## 
## $`xenobiotic metabolism (M164)`
## [1] 10
## 
## $`enriched in activated dendritic cells (II) (M165)`
## [1] 35
## 
## $`TBA (M166)`
## [1] 17
## 
## $`enriched in cell cycle (M167)`
## [1] 15
## 
## $`enriched in dendritic cells (M168)`
## [1] 19
## 
## $`mitosis (TF motif CCAATNNSNNNGCG) (M169)`
## [1] 16
## 
## $`TBA (M170)`
## [1] 12
## 
## $`heme biosynthesis (I) (M171)`
## [1] 11
## 
## $`enriched for TF motif TTCNRGNNNNTTC (M172)`
## [1] 10
## 
## $`erythrocyte differentiation (M173)`
## [1] 14
## 
## $`TBA (M174)`
## [1] 24
## 
## $`cell development (M175)`
## [1] 13
## 
## $`TBA (M176)`
## [1] 10
## 
## $`TBA (M177.0)`
## [1] 34
## 
## $`TBA (M177.1)`
## [1] 15
## 
## $`enriched for promoter motif NATCACGTGAY (putative SREBF1 targets) (M178)`
## [1] 10
## 
## $`enriched for TF motif PAX3 (M179)`
## [1] 10
## 
## $`TBA (M180)`
## [1] 11
## 
## $`nucleotide metabolism (M181)`
## [1] 10
## 
## $`enriched in DNA interacting proteins (M182)`
## [1] 16
## 
## $`TBA (M183)`
## [1] 10
## 
## $`TBA (M184.0)`
## [1] 18
## 
## $`TBA (M184.1)`
## [1] 11
## 
## $`TBA (M185)`
## [1] 14
## 
## $`TBA (M186)`
## [1] 11
## 
## $`metabolism in mitochondria, peroxisome (M187)`
## [1] 12
## 
## $`TBA (M188)`
## [1] 10
## 
## $`extracellular region cluster (GO) (M189)`
## [1] 16
## 
## $`TBA (M190)`
## [1] 11
## 
## $`transmembrane transport (II) (M191)`
## [1] 18
## 
## $`TBA (M192)`
## [1] 12
## 
## $`TBA (M193)`
## [1] 11
## 
## $`TBA (M194)`
## [1] 14
## 
## $`muscle contraction, SRF targets (M195)`
## [1] 12
## 
## $`platelet activation - actin binding (M196)`
## [1] 17
## 
## $`TBA (M197)`
## [1] 14
## 
## $`TBA (M198)`
## [1] 12
## 
## $`platelet activation & blood coagulation (M199)`
## [1] 13
## 
## $`antigen processing and presentation (M200)`
## [1] 11
## 
## $`TBA (M201)`
## [1] 18
## 
## $`enriched in extracellular matrix & associated proteins (M202)`
## [1] 12
## 
## $`TBA (M203)`
## [1] 13
## 
## $`chaperonin mediated protein folding (I) (M204.0)`
## [1] 14
## 
## $`chaperonin mediated protein folding (II) (M204.1)`
## [1] 10
## 
## $`TBA (M205)`
## [1] 11
## 
## $`Wnt signaling pathway (M206)`
## [1] 14
## 
## $`TBA (M207)`
## [1] 10
## 
## $`TBA (M208)`
## [1] 10
## 
## $`lysosome (M209)`
## [1] 10
## 
## $`extracellular matrix, collagen (M210)`
## [1] 32
## 
## $`TBA (M211)`
## [1] 11
## 
## $`purine nucleotide biosynthesis (M212)`
## [1] 12
## 
## $`regulation of transcription, transcription factors (M213)`
## [1] 21
## 
## $`TBA (M214)`
## [1] 12
## 
## $`small GTPase mediated signal transduction (M215)`
## [1] 16
## 
## $`respiratory electron transport chain (mitochondrion) (M216)`
## [1] 12
## 
## $`TBA (source: B cells) (M217)`
## [1] 10
## 
## $`TBA (M218)`
## [1] 15
## 
## $`respiratory electron transport chain (mitochondrion) (M219)`
## [1] 18
## 
## $`TBA (M220)`
## [1] 10
## 
## $`TBA (M221)`
## [1] 16
## 
## $`heme biosynthesis (II) (M222)`
## [1] 13
## 
## $`enriched in T cells (II) (M223)`
## [1] 13
## 
## $`transmembrane and ion transporters (II) (M224)`
## [1] 11
## 
## $`metabolism of steroids (M225)`
## [1] 10
## 
## $`proteasome (M226)`
## [1] 12
## 
## $`translation initiation (M227)`
## [1] 10
## 
## $`olfactory receptors (M228)`
## [1] 13
## 
## $`TBA (M229)`
## [1] 11
## 
## $`cell cycle, mitotic phase (M230)`
## [1] 12
## 
## $`respiratory electron transport chain (mitochondrion) (M231)`
## [1] 11
## 
## $`enriched for TF motif TNCATNTCCYR (M232)`
## [1] 13
## 
## $`TBA (M233)`
## [1] 15
## 
## $`transcription elongation, RNA polymerase II (M234)`
## [1] 14
## 
## $`mitochondrial cluster (M235)`
## [1] 19
## 
## $`TBA (M236)`
## [1] 14
## 
## $`golgi membrane (II) (M237)`
## [1] 17
## 
## $`respiratory electron transport chain (mitochondrion) (M238)`
## [1] 17
## 
## $`enriched in calcium signaling proteins (M239)`
## [1] 15
## 
## $`chromosome Y linked (M240)`
## [1] 17
## 
## $`TBA (M241)`
## [1] 10
## 
## $`TBA (M242)`
## [1] 11
## 
## $`TBA (M243)`
## [1] 12
## 
## $`TBA (M244)`
## [1] 13
## 
## $`translation initiation factor 3 complex (M245)`
## [1] 13
## 
## $`TBA (M246)`
## [1] 16
## 
## $`enriched in nuclear pore complex interacting proteins (M247)`
## [1] 14
## 
## $`TBA (M248)`
## [1] 11
## 
## $`TBA (M249)`
## [1] 10
## 
## $`spliceosome (M250)`
## [1] 12
## 
## $`T cell surface signature (S0)`
## [1] 27
## 
## $`NK cell surface signature (S1)`
## [1] 48
## 
## $`B cell surface signature (S2)`
## [1] 169
## 
## $`Plasma cell surface signature (S3)`
## [1] 24
## 
## $`Monocyte surface signature (S4)`
## [1] 94
## 
## $`DC surface signature (S5)`
## [1] 82
## 
## $`CD4 T cell surface signature Th1-stimulated (S6)`
## [1] 9
## 
## $`CD4 T cell surface signature Th2-stimulated (S7)`
## [1] 13
## 
## $`Naive B cell surface signature (S8)`
## [1] 54
## 
## $`Memory B cell surface signature (S9)`
## [1] 37
## 
## $`Resting dendritic cell surface signature (S10)`
## [1] 75
## 
## $`Activated (LPS) dendritic cell surface signature (S11)`
## [1] 37

# Get all unique genes
all_genes <- sort(unique(unlist(btm)))

# Create the binary matrix
binary_matrix_btm <- sapply(all_genes, function(g) {
    as.numeric(sapply(btm, function(x) g %in% x))
})

# Transpose and convert to a data frame for better visualization
#binary_matrix_btm <- as.data.frame(t(binary_matrix_btm))
binary_matrix_btm[1:5, 1:5]

##      A1CF A2BP1 A2M ABCA1 ABCA13
## [1,]    0     0   0     0      0
## [2,]    0     0   0     0      0
## [3,]    0     0   0     0      0
## [4,]    0     0   0     0      0
## [5,]    0     0   0     0      0

# Add row names
rownames(binary_matrix_btm) <- names(btm)

# Display the binary matrix
binary_matrix_btm[1:5, 1:10]

##                                                A1CF A2BP1 A2M ABCA1 ABCA13
## targets of FOSL1/2 (M0)                           0     0   0     0      0
## integrin cell surface interactions (I) (M1.0)     0     0   0     0      0
## integrin cell surface interactions (II) (M1.1)    0     0   0     0      0
## extracellular matrix (I) (M2.0)                   0     0   0     0      0
## extracellular matrix (II) (M2.1)                  0     0   0     0      0
##                                                ABCA6 ABCA8 ABCB4 ABCB6 ABCC13
## targets of FOSL1/2 (M0)                            0     0     0     0      0
## integrin cell surface interactions (I) (M1.0)      0     0     0     0      0
## integrin cell surface interactions (II) (M1.1)     0     0     0     0      0
## extracellular matrix (I) (M2.0)                    0     0     0     0      0
## extracellular matrix (II) (M2.1)                   0     0     0     0      0

dim(binary_matrix)

## [1]   50 4384

enrichment4.parametric.pos <- run_enrichment(MOFAobject3,
                                             view = "geneExp", factors = 11,
                                             feature.sets = as.matrix(binary_matrix_btm),
                                             sign = "positive",
                                             statistical.test = "parametric"
)

## Intersecting features names in the model and the gene set annotation results in a total of 1277 features.

## 
## Running feature set Enrichment Analysis with the following options...
## View: geneExp 
## Number of feature sets: 124 
## Set statistic: mean.diff 
## Statistical test: parametric

## Subsetting weights with positive sign

## 

plot_enrichment(enrichment4.parametric.pos, 
                factor = 'Factor11', 
                max.pathways = 15
)

plot_enrichment_detailed(enrichment4.parametric.pos, 
                         factor = 'Factor11', 
                         max.genes = 8, 
                         max.pathways = 5
)

## Warning: ggrepel: 7 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

enrichment4.parametric.neg <- run_enrichment(MOFAobject3,
                                             view = "geneExp", factors = 11,
                                             feature.sets = as.matrix(binary_matrix_btm),
                                             sign = "negative",
                                             statistical.test = "parametric"
)

## Intersecting features names in the model and the gene set annotation results in a total of 1277 features.

## 
## Running feature set Enrichment Analysis with the following options...
## View: geneExp 
## Number of feature sets: 124 
## Set statistic: mean.diff 
## Statistical test: parametric

## Subsetting weights with negative sign

## 

plot_enrichment(enrichment4.parametric.neg, 
                factor = 'Factor11', 
                max.pathways = 15
)

plot_enrichment_detailed(enrichment4.parametric.neg, 
                         factor = 'Factor11', 
                         max.genes = 8, 
                         max.pathways = 5
)

## Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

other modules

enrichment.parametric <- run_enrichment(MOFAobject3,
                                             view = "geneExp", factors = c(1:15),
                                             feature.sets = as.matrix(binary_matrix_btm),
                                             sign = "positive",
                                             statistical.test = "parametric"
)

## Intersecting features names in the model and the gene set annotation results in a total of 1277 features.

## 
## Running feature set Enrichment Analysis with the following options...
## View: geneExp 
## Number of feature sets: 124 
## Set statistic: mean.diff 
## Statistical test: parametric

## Subsetting weights with positive sign

## 

weights <- get_weights(MOFAobject3, views = "all", factors = "Factor11", scale = F)
weights.scale  <- get_weights(MOFAobject3, views = "geneExp", factors = "Factor11", scale = T)

a= weights$ab
ggplot(a, aes(x=Factor11)) +
    geom_density() 

a['IgG_PT', ]

## [1] 0.01339264

#weights.scale$geneExp[c('CCL3', 'CXCL8', 'IL1B', 'CCL4'), ]

a = a[order(a[,1]), ]
head(a)

##     IgG1_OVA     IgG3_OVA     IgG1_PRN      IgG1_TT      IgG1_DT     IgG4_OVA 
## -0.003788381 -0.001377719  0.003559370  0.003716228  0.004172358  0.004277007

tail(a)

## IgG2_FIM2/3     IgG4_TT    IgG4_FHA IgG4_FIM2/3     IgG4_PT    IgG4_PRN 
##  0.05131217  0.05355578  0.05548806  0.07388952  0.09171806  0.11257229

d = as.data.frame(tail(names(a), n=200))
#write_tsv(d, file='data/CCL3_MOFA_F3.tsv', col_names = T)
d

##    tail(names(a), n = 200)
## 1                 IgG1_OVA
## 2                 IgG3_OVA
## 3                 IgG1_PRN
## 4                  IgG1_TT
## 5                  IgG1_DT
## 6                 IgG4_OVA
## 7                  IgG_PRN
## 8                 IgG2_OVA
## 9                 IgG1_FHA
## 10                 IgG2_DT
## 11                 IgG_FHA
## 12                 IgG3_DT
## 13                  IgG_PT
## 14                 IgG4_DT
## 15             IgG1_FIM2/3
## 16                 IgG3_PT
## 17                 IgG3_TT
## 18                IgG2_FHA
## 19             IgG3_FIM2/3
## 20                IgG2_PRN
## 21                 IgG2_PT
## 22                 IgG1_PT
## 23                IgG3_FHA
## 24                IgG3_PRN
## 25                 IgG2_TT
## 26             IgG2_FIM2/3
## 27                 IgG4_TT
## 28                IgG4_FHA
## 29             IgG4_FIM2/3
## 30                 IgG4_PT
## 31                IgG4_PRN

cytokineLs

plot_top_weights(MOFAobject3,
                 view = "cytokineL",
                 factor = 11,
                 nfeatures = 5
)

plot_top_weights(MOFAobject3,
                 view = "cell_freq",
                 factor = 11,
                 nfeatures = 5
)

plot_data_scatter(MOFAobject,
                  view = "ab",         # view of interest
                  factor = 11,             # factor of interest
                  features = 11,           # number of features to plot (they are selected by weight)
                  add_lm = TRUE,         # add linear regression
                  color_by = "Meta.infancy_vac"
)

end

Sys.Date()

## [1] "2024-12-13"

getwd()

## [1] "/Users/nthrupp/Documents/Projects/CMI-PB3/models/MOFA_4assays_split2/MOFA"

sessionInfo()

## R version 4.3.1 (2023-06-16)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Sonoma 14.4.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Los_Angeles
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] MOFAdata_1.18.0  MOFA2_1.12.1     lubridate_1.9.3  forcats_1.0.0   
##  [5] stringr_1.5.1    dplyr_1.1.4      purrr_1.0.2      readr_2.1.5     
##  [9] tidyr_1.3.1      tibble_3.2.1     ggplot2_3.5.1    tidyverse_2.0.0 
## [13] BiocStyle_2.30.0
## 
## loaded via a namespace (and not attached):
##   [1] rlang_1.1.4           magrittr_2.0.3        matrixStats_1.4.1    
##   [4] compiler_4.3.1        mgcv_1.9-1            dir.expiry_1.10.0    
##   [7] png_0.1-8             vctrs_0.6.5           reshape2_1.4.4       
##  [10] pkgconfig_2.0.3       crayon_1.5.3          fastmap_1.2.0        
##  [13] backports_1.5.0       magick_2.8.4          XVector_0.42.0       
##  [16] labeling_0.4.3        utf8_1.2.4            rmarkdown_2.27       
##  [19] tzdb_0.4.0            bit_4.0.5             tinytex_0.52         
##  [22] xfun_0.47             zlibbioc_1.48.2       cachem_1.1.0         
##  [25] jsonlite_1.8.8        highr_0.11            rhdf5filters_1.14.1  
##  [28] DelayedArray_0.28.0   BiocParallel_1.36.0   Rhdf5lib_1.24.2      
##  [31] broom_1.0.6           parallel_4.3.1        R6_2.5.1             
##  [34] bslib_0.8.0           stringi_1.8.4         RColorBrewer_1.1-3   
##  [37] reticulate_1.39.0     car_3.1-2             jquerylib_0.1.4      
##  [40] Rcpp_1.0.13           bookdown_0.40         knitr_1.48           
##  [43] IRanges_2.36.0        Matrix_1.6-5          splines_4.3.1        
##  [46] timechange_0.3.0      tidyselect_1.2.1      rstudioapi_0.16.0    
##  [49] abind_1.4-8           yaml_2.3.10           codetools_0.2-20     
##  [52] lattice_0.22-6        plyr_1.8.9            basilisk.utils_1.14.1
##  [55] withr_3.0.1           evaluate_1.0.0        Rtsne_0.17           
##  [58] pillar_1.9.0          BiocManager_1.30.23   ggpubr_0.6.0         
##  [61] filelock_1.0.3        MatrixGenerics_1.14.0 carData_3.0-5        
##  [64] corrplot_0.92         stats4_4.3.1          generics_0.1.3       
##  [67] vroom_1.6.5           S4Vectors_0.40.2      hms_1.1.3            
##  [70] munsell_0.5.1         scales_1.3.0          glue_1.7.0           
##  [73] pheatmap_1.0.12       tools_4.3.1           data.table_1.16.0    
##  [76] fgsea_1.28.0          ggsignif_0.6.4        fastmatch_1.1-4      
##  [79] cowplot_1.1.3         rhdf5_2.46.1          grid_4.3.1           
##  [82] colorspace_2.1-1      nlme_3.1-166          basilisk_1.14.3      
##  [85] HDF5Array_1.30.1      cli_3.6.3             fansi_1.0.6          
##  [88] S4Arrays_1.2.1        uwot_0.2.2            gtable_0.3.5         
##  [91] rstatix_0.7.2         sass_0.4.9            digest_0.6.36        
##  [94] BiocGenerics_0.48.1   SparseArray_1.2.4     ggrepel_0.9.6        
##  [97] farver_2.1.2          htmltools_0.5.8.1     lifecycle_1.0.4      
## [100] bit64_4.0.5